---
title: "MTPS_Bivariate"
author: "Edmund Kelly"
date: "2023-04-30"
output: html_document
---
```{r, setup}
library(rlang)
library(OpenMx)
library(umx)
library(haven)
library(tidyverse)
library(stargazer)
library(eeptools)
library(sjPlot)
library(ggplot2)
library(broom)
library(RColorBrewer)
library(psych)
set.seed(1234)

knitr::opts_knit$set(root.dir = "C:/Users/edmun/Documents/2023-2026 - DPhil Politics/2 - OTHER PAPERS/2 - UNDER REVIEW/UR3 - Twins and trust/Replication materials/Minnesota")
setwd("C:/Users/edmun/Documents/2023-2026 - DPhil Politics/2 - OTHER PAPERS/2 - UNDER REVIEW/UR3 - Twins and trust/Replication materials/Minnesota")

```

```{r Packages, echo = FALSE} 
library(OpenMx)
library(umx) 
library(haven)
library(tidyverse)
library(stargazer)
library(eeptools)
library(sjPlot)
library(ggplot2)
library(broom)
library(RColorBrewer)
library(psych)
library(ggh4x)
library(ggpubr)
set.seed(1234)
```

```{r Pre-Processing}

#Loading data
MTPS <- data.frame(read_spss("MTPS.sav"))


#Removing missing values and incomplete pairs
MTPS[MTPS == 999] <- NA
MTPS[MTPS == 99999999] <- NA
MTPS <- subset(MTPS, complete_pair == 1)

#Converting variables to factors
#MTPS[ , c(6:263, 268:271)] <- lapply(MTPS[ , c(6:263, 268:271)], as.factor)
MTPS$rsex <- as.numeric(MTPS$rsex)
MTPS$BDATE <- as.Date(MTPS$BDATE)
a = as.Date('2024-10-04')
MTPS$AGE <- floor(age_calc(as.Date(MTPS$BDATE), enddate = a, units = "years", precise = TRUE)) - 14

#Turnout

MTPS$turnout04 <- as.factor(MTPS$Qq39)
MTPS <- MTPS %>% mutate(turnout04 = 
                          scale(dplyr::recode(turnout04, 
                                 "1" = 1,
                                 "2" = 0, 
                                 "3" = 0)))

MTPS$lifeturnout <- as.numeric(scale(6-MTPS$Qq41))

MTPS$overallturnout = as.numeric(scale(MTPS$turnout04 + MTPS$lifeturnout))

#Behavioral engagement

MTPS$rally <- as.factor(MTPS$Qq43_A_1)
MTPS <- MTPS %>% mutate(rally = 
                          as.numeric(scale(dplyr::recode(rally, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$workpol <- as.factor(MTPS$Qq43_A_2)
MTPS <- MTPS %>% mutate(workpol = 
                          as.numeric(scale(dplyr::recode(workpol, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$donate <- as.factor(MTPS$Qq43_A_3)
MTPS <- MTPS %>% mutate(donate = 
                          as.numeric(scale(dplyr::recode(donate, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$writepol <- as.factor(MTPS$Qq43_A_5)
MTPS <- MTPS %>% mutate(writepol = 
                          as.numeric(scale(dplyr::recode(writepol, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$behavioral <- as.numeric(scale(MTPS$donate + MTPS$writepol + MTPS$workpol + MTPS$rally))

#Political trust

MTPS$Qq27 = as.numeric(MTPS$Qq27)

MTPS$trustgov <- ((3-MTPS$Qq28))/2

MTPS$biginterest <- (MTPS$Qq29-1)

MTPS$trustpol <- (MTPS$Qq31-1)/2

MTPS$wastetax <- (MTPS$Qq30-1)/2

MTPS$trustgov_scaled = scale(MTPS$trustgov)

MTPS$biginterest_scaled = scale(MTPS$biginterest)

MTPS$wastetax_scaled = scale(MTPS$wastetax)

MTPS$trustpol_scaled = scale(MTPS$trustpol)

MTPS$trust <- as.numeric(MTPS$trustgov + MTPS$biginterest + MTPS$trustpol + MTPS$wastetax)/4
MTPS$trust_scaled = as.numeric(scale(MTPS$trust))

#Equal environment variables

MTPS$bedroom = as.numeric(MTPS$Qq73_A_1)
MTPS$dressalike = as.numeric(MTPS$Qq73_A_3)
MTPS$sameclass = as.numeric(MTPS$Qq73_A_4)
MTPS$samefriends = as.numeric(MTPS$Qq73_A_2)

#Creating twin ID variable
MTPS$TWINID = MTPS$COMPID %% 100
MTPS$TWINID = as.factor(MTPS$TWINID)
MTPS <- MTPS %>% mutate(TWINID = 
                          dplyr::recode(TWINID, 
                                 "0" = 1,
                                 "1" = 2))

#Stacking data by twin pairs
MTPS_desc = MTPS
MTPS <- umx_long2wide(MTPS, famID = "FAMID", twinID = "TWINID", zygosity = "ZYG")

#Creating MZ and DZ datasets
MTPS_MZ <- subset(MTPS, ZYG == 1)
MTPS_DZ <- subset(MTPS, ZYG == 2)

```

```{r Bivariate Cholesky behavioral scale} 

CHO1 <- umxACE(selDVs = c("trust", "behavioral"), selCovs = c("AGE", "rsex"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, tryHard = "yes", optimizer = "SLSQP")

genCorAlg <- mxAlgebra(cov2cor(top.a%*%t(top.a)), name="geneticCorrelation")
genCorCI <- mxCI("geneticCorrelation")

comCorAlg <- mxAlgebra(cov2cor(top.c%*%t(top.c)), name="comCorrelation")
comCorCI <- mxCI("comCorrelation")

envCorAlg <- mxAlgebra(cov2cor(top.e%*%t(top.e)), name="envCorrelation")
envCorCI <- mxCI("envCorrelation")

CHO1 <- mxModel(CHO1, genCorAlg, genCorCI, envCorAlg, envCorCI, comCorAlg, comCorCI)

CHO1 <- mxRun(CHO1, intervals = T)

summary(CHO1, verbose = T)
```

```{r Bivariate Cholesky turnout scale} 

CHO1 <- umxACE(selDVs = c("trust", "overallturnout"), selCovs = c("AGE", "rsex"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, tryHard = "yes", optimizer = "SLSQP")

genCorAlg <- mxAlgebra(cov2cor(top.a%*%t(top.a)), name="geneticCorrelation")
genCorCI <- mxCI("geneticCorrelation")

comCorAlg <- mxAlgebra(cov2cor(top.c%*%t(top.c)), name="comCorrelation")
comCorCI <- mxCI("comCorrelation")

envCorAlg <- mxAlgebra(cov2cor(top.e%*%t(top.e)), name="envCorrelation")
envCorCI <- mxCI("envCorrelation")

CHO1 <- mxModel(CHO1, genCorAlg, genCorCI, envCorAlg, envCorCI, comCorAlg, comCorCI)

CHO1 <- mxRun(CHO1, intervals = T)

summary(CHO1, verbose = T)
```


